# replication file
# Nostalgic Deprivation and Populism: Evidence from 19 European Countries 
# Jeremy Ferwerda, Justin Gest, Tyler Reny
# EJPR
# For questions with code, contact Tyler Reny
# ttreny@gmail.com

rm(list=ls())

# read in packages
pacman::p_load(
  ggeffects, haven, tidyverse, tidytext, estimatr, 
  MASS, readxl, sandwich, survey, srvyr,
  simcf, patchwork,wCorr, weights
)

# set working dir
#setwd()

# read in data
df_t <- read_csv("dat_ejpr_ferwerda_gest_reny.csv")

# Figure A1
hist_populism = ggplot(data=df_t, aes(populism_scale)) +
  geom_histogram(binwidth = 0.1, fill='white',color='black') +
  theme_bw() +
  labs(x='Populism Scale', y='')
ggsave(hist_populism, width=3, height=3,
       file='figs/figureA1.pdf')

# Figure A2
hist_dep = bind_rows(tibble(y=df_t$total_dep_s, var='Deprivation Scale'),
                     tibble(y=df_t$social_dep, var='Social Deprivation'),
                     tibble(y=df_t$econ_dep, var='Economic Deprivation'),
                     tibble(y=df_t$power_dep, var='Power Deprivation')) %>%
  ggplot(aes(x=y)) +
  geom_histogram(binwidth = 0.1, fill='white',color='black') +
  theme_bw() +
  labs(x='Deprivation Scale', y='') +
  facet_wrap(~var, ncol=4)
ggsave(hist_dep, width=10, height=2.5,file='figs/figureA2.pdf')

# Figure C1: Correlations by Country
p <- df_t %>%
  group_by(country) %>%
  summarise(cor_soc_pow=wtd.cor(x=social_dep,y=power_dep,weight=weight)[1],
            cor_soc_econ=wtd.cor(x=social_dep,y=econ_dep,weight=weight)[1],
            cor_pow_econ=wtd.cor(x=power_dep,y=econ_dep,weight=weight)[1]) %>%
  pivot_longer(2:4) %>%
  mutate(name=str_replace(name,"cor_soc_pow","Correlation: Social and Power")) %>%
  mutate(name=str_replace(name,"cor_pow_econ","Correlation: Power and Economic ")) %>%
  mutate(name=str_replace(name,"cor_soc_econ","Correlation: Social and Economic ")) %>%
  ggplot(aes(x=reorder(country,value), y=value)) +
  geom_point(shape=21, fill='white')+
  coord_flip() +
  theme_bw() +
  scale_shape_manual(values = c(20,13,23), name='')+ 
  facet_wrap(~name, ncol=3) + xlab("") + ylab("Correlation Coefficient")
p
ggsave(p, width=9, height=5, file='figs/figureC1.pdf')

# Figure 1 Panel A
d1a <- df_t %>%
  group_by(country) %>%
  summarise(y = weighted.mean(soc_loss, na.rm=T, w=weight),
            l = y - 1.96 * sd(soc_loss, na.rm=T)/sqrt(n()),
            u = y + 1.96 * sd(soc_loss, na.rm=T)/sqrt(n())) %>%
  mutate(name = "Social Deprivation")
d2a <- df_t %>% 
  group_by(country) %>%
  summarise(y = weighted.mean(power_loss, na.rm=T, w=weight),
            l = y - 1.96 * sd(power_loss, na.rm=T)/sqrt(n()),
            u = y + 1.96 * sd(power_loss, na.rm=T)/sqrt(n()))%>%
  mutate(name = "Power Deprivation")
d3a <- df_t %>% 
  group_by(country) %>%
  summarise(y = weighted.mean(econ_loss, na.rm=T, w=weight),
            l = y - 1.96 * sd(econ_loss, na.rm=T)/sqrt(n()),
            u = y + 1.96 * sd(econ_loss, na.rm=T)/sqrt(n())) %>%
  mutate(name = 'Economic Deprivation')
plot1 <- bind_rows(d1a, d2a, d3a) %>%
  ggplot(aes(x=reorder(country,y), y=y, ymin=l, ymax=u, shape=name)) +
  geom_errorbar(width=0, position=position_dodge(width=0.4), color='grey') +
  geom_point(position=position_dodge(width=0.4)) +
  coord_flip() +
  theme_bw() +
  labs(x='', y='Proportion indicating Deprivation') +
  scale_shape_manual(values = c(20,13,23), name='') +
  scale_y_continuous(limits=c(0.15,0.65))
plot1

# Figure 1 Panel B

# this code creates estimates for Figure 1 Panel B; cached in boots_sample_descrips2.csv
# set.seed(12345)
# boots <- rsample::bootstraps(df_t, times = 10000, apparent = TRUE)
# 
# diffInMeans <- function(df_t){
#   full <- df_t %>%
#     summarise(female = weighted.mean(female, na.rm=T, w=weighttot),
#               professionals = weighted.mean(professionals, na.rm=T, w=weighttot),
#               age_18_34 = weighted.mean(age_cut == '18 - 34 yo', na.rm=T, w=weighttot),
#               age_35_49 = weighted.mean(age_cut == '35 - 49 yo', na.rm=T, w=weighttot),
#               age_50_59 = weighted.mean(age_cut == '50 - 59 yo', na.rm=T, w=weighttot),
#               age_60_more = weighted.mean(age_cut == '60 yo and more', na.rm=T, w=weighttot),
#               own_home = weighted.mean(own_home, na.rm=T, w=weighttot),
#               edu_low = weighted.mean(edu_low, na.rm=T, w=weighttot),
#               edu_med = weighted.mean(edu_med, na.rm=T, w=weighttot),
#               edu_high = weighted.mean(edu_high, na.rm=T, w=weighttot),
#               inc = weighted.mean(inc, na.rm=T, w=weighttot),
#               married = weighted.mean(married, na.rm=T, w=weighttot),
#               religiosity = weighted.mean(religiosity, na.rm=T, w=weighttot),
#               rural = weighted.mean(rural, na.rm=T, w=weighttot),
#               foreign_born = weighted.mean(foreign_born, na.rm=T, w=weighttot))%>%
#     pivot_longer(1:15) %>%
#     mutate(var='Full Sample')
#   
#   social_dep <- df_t %>%
#     filter(social_dep >= social_dep_quant) %>%
#     summarise(female = weighted.mean(female, na.rm=T, w=weighttot),
#               professionals = weighted.mean(professionals, na.rm=T, w=weighttot),
#               age_18_34 = weighted.mean(age_cut == '18 - 34 yo', na.rm=T, w=weighttot),
#               age_35_49 = weighted.mean(age_cut == '35 - 49 yo', na.rm=T, w=weighttot),
#               age_50_59 = weighted.mean(age_cut == '50 - 59 yo', na.rm=T, w=weighttot),
#               age_60_more = weighted.mean(age_cut == '60 yo and more', na.rm=T, w=weighttot),
#               own_home = weighted.mean(own_home, na.rm=T, w=weighttot),
#               edu_low = weighted.mean(edu_low, na.rm=T, w=weighttot),
#               edu_med = weighted.mean(edu_med, na.rm=T, w=weighttot),
#               edu_high = weighted.mean(edu_high, na.rm=T, w=weighttot),
#               inc = weighted.mean(inc, na.rm=T, w=weighttot),
#               married = weighted.mean(married, na.rm=T, w=weighttot),
#               religiosity = weighted.mean(religiosity, na.rm=T, w=weighttot),
#               rural = weighted.mean(rural, na.rm=T, w=weighttot),
#               foreign_born = weighted.mean(foreign_born, na.rm=T, w=weighttot)) %>%
#     pivot_longer(1:15) %>%
#     mutate(var='Social Deprivation')
#   
#   
#   econ_dep <- df_t %>%
#     filter(econ_dep >= econ_dep_quant) %>%
#     summarise(female = weighted.mean(female, na.rm=T, w=weighttot),
#               professionals = weighted.mean(professionals, na.rm=T, w=weighttot),
#               age_18_34 = weighted.mean(age_cut == '18 - 34 yo', na.rm=T, w=weighttot),
#               age_35_49 = weighted.mean(age_cut == '35 - 49 yo', na.rm=T, w=weighttot),
#               age_50_59 = weighted.mean(age_cut == '50 - 59 yo', na.rm=T, w=weighttot),
#               age_60_more = weighted.mean(age_cut == '60 yo and more', na.rm=T, w=weighttot),
#               own_home = weighted.mean(own_home, na.rm=T, w=weighttot),
#               edu_low = weighted.mean(edu_low, na.rm=T, w=weighttot),
#               edu_med = weighted.mean(edu_med, na.rm=T, w=weighttot),
#               edu_high = weighted.mean(edu_high, na.rm=T, w=weighttot),
#               inc = weighted.mean(inc, na.rm=T, w=weighttot),
#               married = weighted.mean(married, na.rm=T, w=weighttot),
#               religiosity = weighted.mean(religiosity, na.rm=T, w=weighttot),
#               rural = weighted.mean(rural, na.rm=T, w=weighttot),
#               foreign_born = weighted.mean(foreign_born, na.rm=T, w=weighttot)) %>%
#     pivot_longer(1:15) %>%
#     mutate(var='Economic Deprivation')
#   
#   power_dep <- df_t %>%
#     filter(power_dep >= power_dep_quant) %>%
#     summarise(female = weighted.mean(female, na.rm=T, w=weighttot),
#               professionals = weighted.mean(professionals, na.rm=T, w=weighttot),
#               age_18_34 = weighted.mean(age_cut == '18 - 34 yo', na.rm=T, w=weighttot),
#               age_35_49 = weighted.mean(age_cut == '35 - 49 yo', na.rm=T, w=weighttot),
#               age_50_59 = weighted.mean(age_cut == '50 - 59 yo', na.rm=T, w=weighttot),
#               age_60_more = weighted.mean(age_cut == '60 yo and more', na.rm=T, w=weighttot),
#               own_home = weighted.mean(own_home, na.rm=T, w=weighttot),
#               edu_low = weighted.mean(edu_low, na.rm=T, w=weighttot),
#               edu_med = weighted.mean(edu_med, na.rm=T, w=weighttot),
#               edu_high = weighted.mean(edu_high, na.rm=T, w=weighttot),
#               inc = weighted.mean(inc, na.rm=T, w=weighttot),
#               married = weighted.mean(married, na.rm=T, w=weighttot),
#               religiosity = weighted.mean(religiosity, na.rm=T, w=weighttot),
#               rural = weighted.mean(rural, na.rm=T, w=weighttot),
#               foreign_born = weighted.mean(foreign_born, na.rm=T, w=weighttot)) %>%
#     pivot_longer(1:15) %>%
#     mutate(var='Power Deprivation')
#   
#   social_dep <- left_join(full, social_dep, by='name') %>%
#     mutate(value = value.y - value.x) %>%
#     select(name, value, var.y)
#   
#   econ_dep <- left_join(full, econ_dep, by='name') %>%
#     mutate(value = value.y - value.x) %>%
#     select(name, value, var.y)
#   
#   power_dep <- left_join(full, power_dep, by='name') %>%
#     mutate(value = value.y - value.x) %>%
#     select(name, value, var.y)
#   
#   dat <- bind_rows(social_dep,econ_dep,power_dep) %>%
#     mutate(name = case_when(
#       name == 'female' ~ "Female",
#       name == 'professionals' ~ "Professionals",
#       name == 'age_18_34' ~ "Age 18-34",
#       name == 'age_35_49' ~ "Age 35-49",
#       name == 'age_50_59' ~ "Age 50-59",
#       name == 'age_60_more' ~ "Age over 60",
#       name == 'own_home' ~ "Homeowner",
#       name == 'edu_low' ~ "Lower Education",
#       name == 'edu_med' ~ "Medium Education",
#       name == 'edu_high' ~ "Higher Education",
#       name == 'inc' ~ 'Income',
#       name == 'married' ~ "Married",
#       name == 'religiosity' ~ "High Religiosity",
#       name == 'rural' ~ "Rural Residence",
#       name == 'foreign_born' ~ "Foreign Born"
#     ))
#   return(dat$value)
# }
# n = 1000
# bootsout <- matrix(NA, nrow=45, ncol=n)
# for(i in 1:n){
#   bootsout[,i] <- diffInMeans(rsample::analysis(boots$splits[[i]]))
#   print(i)
# }
# 
# dat <- data.frame(
#   mean = apply(bootsout, 1, mean),
#   l = apply(bootsout, 1, function(x) quantile(x, 0.025)),
#   u = apply(bootsout, 1, function(x) quantile(x, 0.975)),
#   name = c("Female", "Professionals", "Age 18-34","Age 35-49", "Age 50-59",
#            "Age over 60","Homeowner","Lower Education","Medium Education",
#            "Higher Education",'Income',"Married", "High Religiosity","Rural Residence","Foreign Born"),
#   var.y = c(
#     rep('Social Deprivation',15),
#     rep('Economic Deprivation',15),
#     rep('Power Deprivation',15)
#   )
# )
#write_csv(dat, file='boots_samp_descrips2.csv')

dat <- read_csv(file='boots_samp_descrips2.csv')
plot2 <-  dat %>% 
  mutate(name = factor(name, level=rev(dat$name[1:15])),
         var.y = case_when(
           var.y == 'Social Deprivation' ~ 'Social',
           var.y == 'Economic Deprivation' ~ 'Economic', 
           var.y == 'Power Deprivation' ~ 'Power' 
         )) %>%
  ggplot(aes(x=name, y=mean, ymin=l, ymax=u, shape=var.y)) +
  geom_hline(yintercept=0, color='red') +
  geom_errorbar(width=0.1, position=position_dodge(width=0.4), color='grey') +
  geom_point(position=position_dodge(width=0.4)) +
  coord_flip() +
  theme_bw() +
  labs(x='', y='Difference from\nFull Sample Mean (ppt)') +
  scale_shape_manual(values = c(20,13,23), name='Deprivation') +
  scale_y_continuous(limits=c(-0.1, 0.05), 
                     breaks = c(-0.10, -0.05, 0, 0.05),
                     labels = c(-10, -5, 0, 5))
plot2

plot1 <- plot1 + theme(legend.position='none')
p <- plot1 + plot2 
p
ggsave('figs/figure1.pdf', width=9, height=5, p)

# MODELS: Voting and Support-------
my_des <- df_t  %>% as_survey_design(ids=country,weights=weighttot)

demo1 <- "total_dep_s*incumbent + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"
demo2 <- "social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"
demo3 <- "econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"
demo4 <- "power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"
demo5 <- "social_dep*incumbent + econ_dep*incumbent + power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"

model1 <- svyglm(formula(paste0("populism_scale ~ ",demo1)), my_des, family=gaussian(link='identity'))
model2 <- svyglm(formula(paste0("populism_scale ~ ",demo2)), my_des, family=gaussian(link='identity'))
model3 <- svyglm(formula(paste0("populism_scale ~ ",demo3)), my_des, family=gaussian(link='identity'))
model4 <- svyglm(formula(paste0("populism_scale ~ ",demo4)), my_des, family=gaussian(link='identity'))
model5 <- svyglm(formula(paste0("voted_populist ~ ",demo1)), my_des, family=binomial(link='logit'))
model6 <- svyglm(formula(paste0("voted_populist ~ ",demo2)), my_des, family=binomial(link='logit'))
model7 <- svyglm(formula(paste0("voted_populist ~ ",demo3)), my_des, family=binomial(link='logit'))
model8 <- svyglm(formula(paste0("voted_populist ~ ",demo4)), my_des, family=binomial(link='logit'))

texreg::texreg(list(model1, model2, model3, model4, model5, model6, model7, model8),                   
               omit.coef="(deu)|(ita)|(esp)|(nld)|(fra)|(swe)|(rou)|(hun)|(pol)|(est)|(aut)|(dnk)|(cze)|(lva)|(ltu)|(svn)|(svk)|(bgr)",
               naive=T,reorder.coef = c(1,2,21,23,25,3:19,20,22,24,26),
               custom.coef.names = c('Intercept','Deprivation Scale','Incumbent','Female','Professional',
                                     'Service Worker','Worker','Age 35-49','Age 50-59','Age 60 plus',
                                     'Med Income','High Income','Med Education','High Education','Married',
                                     'Foreign Born','Worried COVID','Unemployed','Religious','Dep Scale * Incumbent',
                                     'Social Deprivation','Soc Dep * Incumbent','Economic Deprivation', 'Econ Dep * Incumbent',
                                     'Power Deprivation','Power Dep * Incumbent'),custom.note = '',
               custom.header = list('Populist Attitudes'=1:4,'Populist Voting'=5:8),
               file = 'tabs/tableB1.tex')

# all dep vars in same model
model9 <- svyglm(formula(paste0("populism_scale ~ ",demo5)), my_des, family=gaussian(link='identity'))
model10 <- svyglm(formula(paste0("voted_populist ~ ",demo5)), my_des, family=binomial(link='logit'))

texreg::texreg(list(model9, model10),                   
               omit.coef="(deu)|(ita)|(esp)|(nld)|(fra)|(swe)|(rou)|(hun)|(pol)|(est)|(aut)|(dnk)|(cze)|(lva)|(ltu)|(svn)|(svk)|(bgr)",
               naive=T,
               custom.coef.names = c('Intercept','Social Deprivation','Incumbent',
                                     'Economic Deprivation','Power Deprivation',
                                     'Female','Professional',
                                     'Service Worker','Worker','Age 35-49','Age 50-59','Age 60 plus',
                                     'Med Income','High Income','Med Education','High Education','Married',
                                     'Foreign Born','Worried COVID','Unemployed','Religious',
                                     'Soc Dep * Incumbent','Econ Dep * Incumbent','Power Dep * Incumbent'),
               reorder.coef = c(1,2,4,5, 3,6:24),
               custom.note = '',
               custom.header = list('Populist Attitudes'=1,'Populist Voting'=2),
               file = 'tabs/tableB3.tex')

# PD predictive of non voting
model11a <- svyglm(formula(paste0("didntvote ~ ",demo4)), my_des, family=binomial(link='logit'))
model11b <- svyglm(formula(paste0("didntvote ~ ",demo1)), my_des, family=binomial(link='logit'))

# vote center left
countries <- c('Austria (AUT)','Denmark (DNK)',
               'France (FRA)', 'Germany (DEU)',
               'Italy (ITA)','Netherlands (NLD)',
               'Spain (ESP)','Sweden (SWE)',
               'United Kingdom (GBR)')
my_des_cl <- df_t  %>% 
  filter(country %in% countries) %>%
  as_survey_design(ids=country,weights=weighttot)

demo_cl <- "econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ deu + ita + fra"
demo_cl2 <- "total_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ deu + ita + fra"

model20a <- svyglm(formula(paste0("voted_center_left_manual ~ ",demo_cl)), my_des_cl, family=binomial(link='logit'))
model20b <- svyglm(formula(paste0("voted_center_left_manual ~ ",demo_cl2)), my_des_cl, family=binomial(link='logit'))

texreg::texreg(list(model11a,model11b, model20a, model20b),                   
               omit.coef="(deu)|(ita)|(esp)|(nld)|(fra)|(swe)|(rou)|(hun)|(pol)|(est)|(aut)|(dnk)|(cze)|(lva)|(ltu)|(svn)|(svk)|(bgr)",
               naive=T,
               custom.coef.names = c('Intercept','Power Deprivation','Incumbent',
                                     'Female','Professional',
                                     'Service Worker','Worker','Age 35-49','Age 50-59','Age 60 plus',
                                     'Med Income','High Income','Med Education','High Education','Married',
                                     'Foreign Born','Worried COVID','Unemployed','Religious',
                                     'Power Dep * Incumbent','Total Dep','Total Dep * Incumbent',
                                     'Econ Deprivation','Econ Dep * Incumbent',
                                     'Total Dep','Total Dep * Incumbent'),
               reorder.coef = c(1,2,21,23,3:20,
                                22,24),
               custom.note = '',
               custom.header = list('Non-Voting'=1:2,'Vote Social Dem'=3:4),
               file = 'tabs/tableB6.tex')

# make plot
set.seed(12345)
pe <- coef(model1)
vc <- vcov(model1)
simbetas1 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model2)
vc <- vcov(model2)
simbetas2 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model3)
vc <- vcov(model3)
simbetas3 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model4)
vc <- vcov(model4)
simbetas4 <- MASS::mvrnorm(10000, pe, vc)

mm1 <- cfMake(populism_scale ~ total_dep_s*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze +
                lva + ltu + svn + svk + bgr, data = df_t, nscen = 1)
mm2 <- cfMake(populism_scale ~ social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze +
                lva + ltu + svn + svk + bgr, data = df_t, nscen = 1)
mm3 <- cfMake(populism_scale ~ econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze +
                lva + ltu + svn + svk + bgr, data = df_t, nscen = 1)
mm4 <- cfMake(populism_scale ~ power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze +
                lva + ltu + svn + svk + bgr, data = df_t, nscen = 1)

mm1 <- cfChange(mm1, 'total_dep_s',xpre=0, x=1, scen=1)
mm1 <- cfChange(mm1, 'incumbent',xpre=0, x=0, scen=1)

mm2 <- cfChange(mm2, 'social_dep',xpre=0, x=1, scen=1)
mm2 <- cfChange(mm2, 'incumbent',xpre=0, x=0, scen=1)

mm3 <- cfChange(mm3, 'econ_dep',xpre=0, x=1, scen=1)
mm3 <- cfChange(mm3, 'incumbent',xpre=0, x=0, scen=1)

mm4 <- cfChange(mm4, 'power_dep',xpre=0, x=1, scen=1)
mm4 <- cfChange(mm4, 'incumbent',xpre=0, x=0, scen=1)


out1 <- linearsimfd(mm1, simbetas1, ci=0.95) %>% as.data.frame() %>% mutate(iv='Deprivation Scale', dv='Populist Attitudes')
out2 <- linearsimfd(mm2, simbetas2, ci=0.95) %>% as.data.frame() %>% mutate(iv='Social Deprivation', dv='Populist Attitudes')
out3 <- linearsimfd(mm3, simbetas3, ci=0.95) %>% as.data.frame() %>% mutate(iv='Economic Deprivation', dv='Populist Attitudes')
out4 <- linearsimfd(mm4, simbetas4, ci=0.95) %>% as.data.frame() %>% mutate(iv='Power Deprivation', dv='Populist Attitudes')
populist_attitudes_out <- bind_rows(out1, out2, out3, out4) 

# voting
set.seed(12345)
pe <- coef(model5)
vc <- vcov(model5)
simbetas5 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model6)
vc <- vcov(model6)
simbetas6 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model7)
vc <- vcov(model7)
simbetas7 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model8)
vc <- vcov(model8)
simbetas8 <- MASS::mvrnorm(10000, pe, vc)

mm1 <- cfMake(voted_populist ~ total_dep_s*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t, nscen = 1)
mm2 <- cfMake(voted_populist ~ social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t, nscen = 1)
mm3 <- cfMake(voted_populist ~ econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t, nscen = 1)
mm4 <- cfMake(voted_populist ~ power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t, nscen = 1)

mm1 <- cfChange(mm1, 'total_dep_s',xpre=0, x=1, scen=1)
mm1 <- cfChange(mm1, 'incumbent',xpre=0, x=0, scen=1)

mm2 <- cfChange(mm2, 'social_dep',xpre=0, x=1, scen=1)
mm2 <- cfChange(mm2, 'incumbent',xpre=0, x=0, scen=1)

mm3 <- cfChange(mm3, 'econ_dep',xpre=0, x=1, scen=1)
mm3 <- cfChange(mm3, 'incumbent',xpre=0, x=0, scen=1)

mm4 <- cfChange(mm4, 'power_dep',xpre=0, x=1, scen=1)
mm4 <- cfChange(mm4, 'incumbent',xpre=0, x=0, scen=1)


out1 <- logitsimfd(mm1, simbetas5, ci=0.95) %>% as.data.frame() %>% mutate(iv='Deprivation Scale', dv='Populist Voting')
out2 <- logitsimfd(mm2, simbetas6, ci=0.95) %>% as.data.frame() %>% mutate(iv='Social Deprivation', dv='Populist Voting')
out3 <- logitsimfd(mm3, simbetas7, ci=0.95) %>% as.data.frame() %>% mutate(iv='Economic Deprivation', dv='Populist Voting')
out4 <- logitsimfd(mm4, simbetas8, ci=0.95) %>% as.data.frame() %>% mutate(iv='Power Deprivation', dv='Populist Voting')
populist_voting_out <- bind_rows(out1, out2, out3, out4) 

# left and right disaggregated
# vote populist left and right manual
countries <- c('Austria (AUT)','Czech Republic (CZE)', 'Denmark (DNK)',
               'France (FRA)', 'Germany (DEU)','Italy (ITA)','Netherlands (NLD)')
my_des_fl <- df_t  %>% 
  filter(country %in% countries) %>%
  as_survey_design(ids=country,weights=weighttot)

model12 <- svyglm(formula(paste0("voted_far_left_manual ~ ",demo1)), my_des_fl, family=binomial(link='logit'))
model13 <- svyglm(formula(paste0("voted_far_left_manual ~ ",demo2)), my_des_fl, family=binomial(link='logit'))
model14 <- svyglm(formula(paste0("voted_far_left_manual ~ ",demo3)), my_des_fl, family=binomial(link='logit'))
model15 <- svyglm(formula(paste0("voted_far_left_manual ~ ",demo4)), my_des_fl, family=binomial(link='logit'))

model16 <- svyglm(formula(paste0("voted_far_right_manual ~ ",demo1)), my_des, family=binomial(link='logit'))
model17 <- svyglm(formula(paste0("voted_far_right_manual ~ ",demo2)), my_des, family=binomial(link='logit'))
model18 <- svyglm(formula(paste0("voted_far_right_manual ~ ",demo3)), my_des, family=binomial(link='logit'))
model19 <- svyglm(formula(paste0("voted_far_right_manual ~ ",demo4)), my_des, family=binomial(link='logit'))
texreg::texreg(list(model12, model13, model14, model15, model16, model17,model18,model19),
               omit.coef="(deu)|(ita)|(esp)|(nld)|(fra)|(swe)|(rou)|(hun)|(pol)|(est)|(aut)|(dnk)|(cze)|(lva)|(ltu)|(svn)|(svk)|(bgr)",
               naive=T,reorder.coef = c(1,2,21,23,25,3:19,20,22,24,26),
               custom.coef.names = c('Intercept','Deprivation Scale','Incumbent','Female','Professional',
                                     'Service Worker','Worker','Age 35-49','Age 50-59','Age 60 plus',
                                     'Med Income','High Income','Med Education','High Education','Married',
                                     'Foreign Born','Worried COVID','Unemployed','Religious','Dep Scale * Incumbent',
                                     'Social Deprivation','Soc Dep * Incumbent','Economic Deprivation', 'Econ Dep * Incumbent',
                                     'Power Deprivation','Power Dep * Incumbent'),custom.note = '',
               custom.header = list('Populist Left Voting'=1:4,'Populist Right Voting'=5:8),
               file = 'tabs/tableB5.tex')

# save PEs for inclusion in Fig 2
# voting
set.seed(12345)
pe <- coef(model12)
vc <- vcov(model12)
simbetas12 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model13)
vc <- vcov(model13)
simbetas13 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model14)
vc <- vcov(model14)
simbetas14 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model15)
vc <- vcov(model15)
simbetas15 <- MASS::mvrnorm(10000, pe, vc)
set.seed(12345)
pe <- coef(model16)
vc <- vcov(model16)
simbetas16 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model17)
vc <- vcov(model17)
simbetas17 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model18)
vc <- vcov(model18)
simbetas18 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model19)
vc <- vcov(model19)
simbetas19 <- MASS::mvrnorm(10000, pe, vc)

mm1 <- cfMake(voted_far_left_manual ~ total_dep_s*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + nld + fra + aut + dnk, data = df_t, nscen = 1)
mm2 <- cfMake(voted_far_left_manual ~ social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + nld + fra + aut + dnk, data = df_t, nscen = 1)
mm3 <- cfMake(voted_far_left_manual ~ econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + nld + fra + aut + dnk, data = df_t, nscen = 1)
mm4 <- cfMake(voted_far_left_manual ~ power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + nld + fra + aut + dnk, data = df_t, nscen = 1)

mm5 <- cfMake(voted_far_right_manual  ~ total_dep_s*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t, nscen = 1)
mm6 <- cfMake(voted_far_right_manual ~ social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t, nscen = 1)
mm7 <- cfMake(voted_far_right_manual ~ econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t, nscen = 1)
mm8 <- cfMake(voted_far_right_manual ~ power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t, nscen = 1)

mm1 <- cfChange(mm1, 'total_dep_s',xpre=0, x=1, scen=1)
mm1 <- cfChange(mm1, 'incumbent',xpre=0, x=0, scen=1)
mm2 <- cfChange(mm2, 'social_dep',xpre=0, x=1, scen=1)
mm2 <- cfChange(mm2, 'incumbent',xpre=0, x=0, scen=1)
mm3 <- cfChange(mm3, 'econ_dep',xpre=0, x=1, scen=1)
mm3 <- cfChange(mm3, 'incumbent',xpre=0, x=0, scen=1)
mm4 <- cfChange(mm4, 'power_dep',xpre=0, x=1, scen=1)
mm4 <- cfChange(mm4, 'incumbent',xpre=0, x=0, scen=1)

mm5 <- cfChange(mm5, 'total_dep_s',xpre=0, x=1, scen=1)
mm5 <- cfChange(mm5, 'incumbent',xpre=0, x=0, scen=1)
mm6 <- cfChange(mm6, 'social_dep',xpre=0, x=1, scen=1)
mm6 <- cfChange(mm6, 'incumbent',xpre=0, x=0, scen=1)
mm7 <- cfChange(mm7, 'econ_dep',xpre=0, x=1, scen=1)
mm7 <- cfChange(mm7, 'incumbent',xpre=0, x=0, scen=1)
mm8 <- cfChange(mm8, 'power_dep',xpre=0, x=1, scen=1)
mm8 <- cfChange(mm8, 'incumbent',xpre=0, x=0, scen=1)

out1 <- logitsimfd(mm1, simbetas12, ci=0.95) %>% as.data.frame() %>% mutate(iv='Deprivation Scale', dv='Populist Voting')
out2 <- logitsimfd(mm2, simbetas13, ci=0.95) %>% as.data.frame() %>% mutate(iv='Social Deprivation', dv='Populist Voting')
out3 <- logitsimfd(mm3, simbetas14, ci=0.95) %>% as.data.frame() %>% mutate(iv='Economic Deprivation', dv='Populist Voting')
out4 <- logitsimfd(mm4, simbetas15, ci=0.95) %>% as.data.frame() %>% mutate(iv='Power Deprivation', dv='Populist Voting')

out5 <- logitsimfd(mm5, simbetas16, ci=0.95) %>% as.data.frame() %>% mutate(iv='Deprivation Scale', dv='Populist Voting')
out6 <- logitsimfd(mm6, simbetas17, ci=0.95) %>% as.data.frame() %>% mutate(iv='Social Deprivation', dv='Populist Voting')
out7 <- logitsimfd(mm7, simbetas18, ci=0.95) %>% as.data.frame() %>% mutate(iv='Economic Deprivation', dv='Populist Voting')
out8 <- logitsimfd(mm8, simbetas19, ci=0.95) %>% as.data.frame() %>% mutate(iv='Power Deprivation', dv='Populist Voting')

populist_voting_out_lr <- bind_rows(out1, out2, out3, out4, out5, out6, out7, out8)  %>%
  mutate(parties = c('Left-wing Parties','Left-wing Parties','Left-wing Parties','Left-wing Parties',
                     'Right-wing Parties','Right-wing Parties','Right-wing Parties','Right-wing Parties'),
         dv = 'Populist Voting')

plot1 <- bind_rows(populist_attitudes_out, populist_voting_out) %>%
  mutate(iv = factor(iv, levels=rev(c('Deprivation Scale','Social Deprivation','Economic Deprivation', 'Power Deprivation')))) %>%
  ggplot(aes(x=iv, y=pe, ymin=lower, ymax=upper)) +
  coord_flip() + 
  geom_hline(yintercept=0, linetype=2, color='red') +
  geom_errorbar(width=0.1, position=position_dodge(width=0.6)) +
  geom_point(size=7, shape=21, fill='white', position=position_dodge(width=0.6)) +
  facet_wrap(~dv) +
  labs(x='',y='Change in DV\n(min-max Deprivation)') +
  theme_bw() +
  geom_text(aes(x=iv, y=pe, label=round(pe,2)),size=2,
            position=position_dodge(width=0.6),
            show.legend = FALSE) 
plot1
ggsave(plot1, width=6, height=3.5, file='figs/figure2.pdf')

# other voting left / right parties
populist_voting_out_lr[1:4,1:3] <- populist_voting_out_lr[1:4,1:3]/sd(df_t$voted_far_left_manual,na.rm=T)
populist_voting_out_lr[5:8,1:3] <- populist_voting_out_lr[5:8,1:3]/sd(df_t$voted_far_right_manual,na.rm=T)

plot2 <- populist_voting_out_lr %>%
  mutate(iv = factor(iv, levels=rev(c('Deprivation Scale','Social Deprivation','Economic Deprivation', 'Power Deprivation')))) %>%
  ggplot(aes(x=iv, y=pe, ymin=lower, ymax=upper, shape=parties, color=parties)) +
  coord_flip() + 
  geom_hline(yintercept=0, linetype=2, color='red') +
  geom_pointrange(position=position_dodge(width=0.6), size=1.5, fill='white') +
  facet_wrap(~dv) +
  labs(x='',y='SD Change in DV\n(min-max Deprivation)') +
  theme_bw() +
  scale_shape_manual(values=c(21,23),name='') +
  scale_color_manual(values=c('blue','red'),name='') +
  geom_text(aes(x=iv, y=pe, label=round(pe,2)),size=2,
            position=position_dodge(width=0.6),
            show.legend = FALSE) 
plot2
ggsave(plot2, width=6, height=3.5, file='figs/figureC2.pdf')

# MODELS ROBUST: Voting and Support-------
demo6 <- "total_dep_s*incumbent + worried_corona + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"
demo7 <- "social_dep*incumbent + worried_corona + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"
demo8 <- "econ_dep*incumbent + worried_corona + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"
demo9 <- "power_dep*incumbent + worried_corona + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"

model1 <- svyglm(formula(paste0("populism_scale ~ ",demo6)), my_des, family=gaussian(link='identity'))
model2 <- svyglm(formula(paste0("populism_scale ~ ",demo7)), my_des, family=gaussian(link='identity'))
model3 <- svyglm(formula(paste0("populism_scale ~ ",demo8)), my_des, family=gaussian(link='identity'))
model4 <- svyglm(formula(paste0("populism_scale ~ ",demo9)), my_des, family=gaussian(link='identity'))
model5 <- svyglm(formula(paste0("voted_populist ~ ",demo6)), my_des, family=binomial(link='logit'))
model6 <- svyglm(formula(paste0("voted_populist ~ ",demo7)), my_des, family=binomial(link='logit'))
model7 <- svyglm(formula(paste0("voted_populist ~ ",demo8)), my_des, family=binomial(link='logit'))
model8 <- svyglm(formula(paste0("voted_populist ~ ",demo9)), my_des, family=binomial(link='logit'))

texreg::texreg(list(model1, model2, model3, model4, model5, model6, model7, model8),
               omit.coef="(deu)|(ita)|(esp)|(nld)|(fra)|(swe)|(rou)|(hun)|(pol)|(est)|(aut)|(dnk)|(cze)|(lva)|(ltu)|(svn)|(svk)|(bgr)",
               naive=T,
               reorder.coef = c(1,2,6,8,10,3,5,7,9,11,4),
               custom.coef.names = c('Intercept','Deprivation Scale','Incumbent',
                                     'Worried COVID','Dep Scale * Incumbent',
                                     'Social Deprivation','Soc Dep * Incumbent',
                                     'Economic Deprivation','Econ Dep * Incumbent',
                                     'Power Deprivation','Power Dep * Incumbent'
               ),custom.note = '',
               custom.header = list('Populist Attitudes'=1:4,'Populist Voting'=5:8),
               file = 'tabs/tableB2.tex')

# MODELS: Left wing vs right wing ----

leftist <- df_t %>% filter(leftist==1) %>% as_survey_design(ids=country,weights=weighttot)
rightist <- df_t %>% filter(rightist==1) %>% as_survey_design(ids=country,weights=weighttot)

model1 <- svyglm(formula(paste0("populism_scale ~ ",demo1)), leftist, family=gaussian(link='identity'))
model2 <- svyglm(formula(paste0("populism_scale ~ ",demo2)), leftist, family=gaussian(link='identity'))
model3 <- svyglm(formula(paste0("populism_scale ~ ",demo3)), leftist, family=gaussian(link='identity'))
model4 <- svyglm(formula(paste0("populism_scale ~ ",demo4)), leftist, family=gaussian(link='identity'))
model5 <- svyglm(formula(paste0("populism_scale ~ ",demo1)), rightist, family=gaussian(link='identity'))
model6 <- svyglm(formula(paste0("populism_scale ~ ",demo2)), rightist, family=gaussian(link='identity'))
model7 <- svyglm(formula(paste0("populism_scale ~ ",demo3)), rightist, family=gaussian(link='identity'))
model8 <- svyglm(formula(paste0("populism_scale ~ ",demo4)), rightist, family=gaussian(link='identity'))
texreg::texreg(list(model1, model2, model3, model4, model5, model6, model7, model8),                   omit.coef="(deu)|(ita)|(esp)|(nld)|(fra)|(swe)|(rou)|(hun)|(pol)|(est)|(aut)|(dnk)|(cze)|(lva)|(ltu)|(svn)|(svk)|(bgr)",
               naive=T,reorder.coef = c(1,2,21,23,25,3:19,20,22,24,26),
               custom.coef.names = c('Intercept','Deprivation Scale','Incumbent','Female','Professional',
                                     'Service Worker','Worker','Age 35-49','Age 50-59','Age 60 plus',
                                     'Med Income','High Income','Med Education','High Education','Married',
                                     'Foreign Born','Worried COVID','Unemployed','Religious','Dep Scale * Incumbent',
                                     'Social Deprivation','Soc Dep * Incumbent','Economic Deprivation', 'Econ Dep * Incumbent',
                                     'Power Deprivation','Power Dep * Incumbent'),custom.note = '',
               custom.header = list('Populist Attitudes (Left-Wing)'=1:4,'Populist Attitudes (Right-Wing)'=5:8),
               file = 'tabs/tableB7.tex')

set.seed(12345)
pe <- coef(model1); vc <- vcov(model1)
simbetas1 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model2); vc <- vcov(model2)
simbetas2 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model3); vc <- vcov(model3)
simbetas3 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model4); vc <- vcov(model4)
simbetas4 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model5); vc <- vcov(model5)
simbetas5 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model6); vc <- vcov(model6)
simbetas6 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model7); vc <- vcov(model7)
simbetas7 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model8); vc <- vcov(model8)
simbetas8 <- MASS::mvrnorm(10000, pe, vc)

mm1 <- cfMake(populism_scale ~ total_dep_s * incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze +lva + ltu + svn + svk + bgr, data = df_t[df_t$leftist==1,], nscen = 1)
mm2 <- cfMake(populism_scale ~ social_dep* incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$leftist==1,], nscen = 1)
mm3 <- cfMake(populism_scale ~ econ_dep* incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$leftist==1,], nscen = 1)
mm4 <- cfMake(populism_scale ~ power_dep* incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$leftist==1,], nscen = 1)
mm5 <- cfMake(populism_scale ~ total_dep_s* incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$rightist==1,], nscen = 1)
mm6 <- cfMake(populism_scale ~ social_dep* incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$rightist==1,], nscen = 1)
mm7 <- cfMake(populism_scale ~ econ_dep* incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$rightist==1,], nscen = 1)
mm8 <- cfMake(populism_scale ~ power_dep* incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$rightist==1,], nscen = 1)

mm1 <- cfChange(mm1, 'total_dep_s',xpre=0, x=1, scen=1)
mm2 <- cfChange(mm2, 'social_dep',xpre=0, x=1, scen=1)
mm3 <- cfChange(mm3, 'econ_dep',xpre=0, x=1, scen=1)
mm4 <- cfChange(mm4, 'power_dep',xpre=0, x=1, scen=1)
mm5 <- cfChange(mm5, 'total_dep_s',xpre=0, x=1, scen=1)
mm6 <- cfChange(mm6, 'social_dep',xpre=0, x=1, scen=1)
mm7 <- cfChange(mm7, 'econ_dep',xpre=0, x=1, scen=1)
mm8 <- cfChange(mm8, 'power_dep',xpre=0, x=1, scen=1)

mm1 <- cfChange(mm1, 'incumbent',xpre=0, x=0, scen=1)
mm2 <- cfChange(mm2, 'incumbent',xpre=0, x=0, scen=1)
mm3 <- cfChange(mm3, 'incumbent',xpre=0, x=0, scen=1)
mm4 <- cfChange(mm4, 'incumbent',xpre=0, x=0, scen=1)
mm5 <- cfChange(mm5, 'incumbent',xpre=0, x=0, scen=1)
mm6 <- cfChange(mm6, 'incumbent',xpre=0, x=0, scen=1)
mm7 <- cfChange(mm7, 'incumbent',xpre=0, x=0, scen=1)
mm8 <- cfChange(mm8, 'incumbent',xpre=0, x=0, scen=1)

out1 <- linearsimfd(mm1, simbetas1, ci=0.95) %>% as.data.frame() %>% mutate(iv='Deprivation Scale', dv='Populist Attitudes', Ideology='Left')
out2 <- linearsimfd(mm2, simbetas2, ci=0.95) %>% as.data.frame() %>% mutate(iv='Social Deprivation', dv='Populist Attitudes', Ideology='Left')
out3 <- linearsimfd(mm3, simbetas3, ci=0.95) %>% as.data.frame() %>% mutate(iv='Economic Deprivation', dv='Populist Attitudes', Ideology='Left')
out4 <- linearsimfd(mm4, simbetas4, ci=0.95) %>% as.data.frame() %>% mutate(iv='Power Deprivation', dv='Populist Attitudes', Ideology='Left')
out5 <- linearsimfd(mm5, simbetas5, ci=0.95) %>% as.data.frame() %>% mutate(iv='Deprivation Scale', dv='Populist Attitudes', Ideology='Right')
out6 <- linearsimfd(mm6, simbetas6, ci=0.95) %>% as.data.frame() %>% mutate(iv='Social Deprivation', dv='Populist Attitudes', Ideology='Right')
out7 <- linearsimfd(mm7, simbetas7, ci=0.95) %>% as.data.frame() %>% mutate(iv='Economic Deprivation', dv='Populist Attitudes', Ideology='Right')
out8 <- linearsimfd(mm8, simbetas8, ci=0.95) %>% as.data.frame() %>% mutate(iv='Power Deprivation', dv='Populist Attitudes', Ideology='Right')

populist_attitudes_out <- bind_rows(out1, out2, out3, out4, out5, out6, out7, out8) 


model1 <- svyglm(formula(paste0("voted_populist ~ ",demo1)), leftist, family=binomial(link='logit'))
model2 <- svyglm(formula(paste0("voted_populist ~ ",demo2)), leftist, family=binomial(link='logit'))
model3 <- svyglm(formula(paste0("voted_populist ~ ",demo3)), leftist, family=binomial(link='logit'))
model4 <- svyglm(formula(paste0("voted_populist ~ ",demo4)), leftist, family=binomial(link='logit'))
model5 <- svyglm(formula(paste0("voted_populist ~ ",demo1)), rightist, family=binomial(link='logit'))
model6 <- svyglm(formula(paste0("voted_populist ~ ",demo2)), rightist, family=binomial(link='logit'))
model7 <- svyglm(formula(paste0("voted_populist ~ ",demo3)), rightist, family=binomial(link='logit'))
model8 <- svyglm(formula(paste0("voted_populist ~ ",demo4)), rightist, family=binomial(link='logit'))
texreg::texreg(list(model1, model2, model3, model4, model5, model6, model7, model8),                   
               omit.coef="(deu)|(ita)|(esp)|(nld)|(fra)|(swe)|(rou)|(hun)|(pol)|(est)|(aut)|(dnk)|(cze)|(lva)|(ltu)|(svn)|(svk)|(bgr)",
               naive=T,reorder.coef = c(1,2,21,23,25,3:19,20,22,24,26),
               custom.coef.names = c('Intercept','Deprivation Scale','Incumbent','Female','Professional',
                                     'Service Worker','Worker','Age 35-49','Age 50-59','Age 60 plus',
                                     'Med Income','High Income','Med Education','High Education','Married',
                                     'Foreign Born','Worried COVID','Unemployed','Religious','Dep Scale * Incumbent',
                                     'Social Deprivation','Soc Dep * Incumbent','Economic Deprivation', 'Econ Dep * Incumbent',
                                     'Power Deprivation','Power Dep * Incumbent'),custom.note = '',
               custom.header = list('Populist Voting (Left-Wing)'=1:4,'Populist Voting (Right-Wing)'=5:8),
               file = 'tabs/tableB8.tex')

set.seed(12345)
pe <- coef(model1); vc <- vcov(model1)
simbetas1 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model2); vc <- vcov(model2)
simbetas2 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model3); vc <- vcov(model3)
simbetas3 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model4); vc <- vcov(model4)
simbetas4 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model5); vc <- vcov(model5)
simbetas5 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model6); vc <- vcov(model6)
simbetas6 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model7); vc <- vcov(model7)
simbetas7 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model8); vc <- vcov(model8)
simbetas8 <- MASS::mvrnorm(10000, pe, vc)

mm1 <- cfMake(voted_populist ~ total_dep_s*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$leftist==1,], nscen = 1)
mm2 <- cfMake(voted_populist ~ social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$leftist==1,], nscen = 1)
mm3 <- cfMake(voted_populist ~ econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$leftist==1,], nscen = 1)
mm4 <- cfMake(voted_populist ~ power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$leftist==1,], nscen = 1)
mm5 <- cfMake(voted_populist ~ total_dep_s*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$rightist==1,], nscen = 1)
mm6 <- cfMake(voted_populist ~ social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$rightist==1,], nscen = 1)
mm7 <- cfMake(voted_populist ~ econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$rightist==1,], nscen = 1)
mm8 <- cfMake(voted_populist ~ power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr, data = df_t[df_t$rightist==1,], nscen = 1)

mm1 <- cfChange(mm1, 'total_dep_s',xpre=0, x=1, scen=1)
mm2 <- cfChange(mm2, 'social_dep',xpre=0, x=1, scen=1)
mm3 <- cfChange(mm3, 'econ_dep',xpre=0, x=1, scen=1)
mm4 <- cfChange(mm4, 'power_dep',xpre=0, x=1, scen=1)
mm5 <- cfChange(mm5, 'total_dep_s',xpre=0, x=1, scen=1)
mm6 <- cfChange(mm6, 'social_dep',xpre=0, x=1, scen=1)
mm7 <- cfChange(mm7, 'econ_dep',xpre=0, x=1, scen=1)
mm8 <- cfChange(mm8, 'power_dep',xpre=0, x=1, scen=1)

mm1 <- cfChange(mm1, 'incumbent',xpre=0, x=0, scen=1)
mm2 <- cfChange(mm2, 'incumbent',xpre=0, x=0, scen=1)
mm3 <- cfChange(mm3, 'incumbent',xpre=0, x=0, scen=1)
mm4 <- cfChange(mm4, 'incumbent',xpre=0, x=0, scen=1)
mm5 <- cfChange(mm5, 'incumbent',xpre=0, x=0, scen=1)
mm6 <- cfChange(mm6, 'incumbent',xpre=0, x=0, scen=1)
mm7 <- cfChange(mm7, 'incumbent',xpre=0, x=0, scen=1)
mm8 <- cfChange(mm8, 'incumbent',xpre=0, x=0, scen=1)

out1 <- logitsimfd(mm1, simbetas1, ci=0.95) %>% as.data.frame() %>% mutate(iv='Deprivation Scale', dv='Populist Voting', Ideology='Left')
out2 <- logitsimfd(mm2, simbetas2, ci=0.95) %>% as.data.frame() %>% mutate(iv='Social Deprivation', dv='Populist Voting', Ideology='Left')
out3 <- logitsimfd(mm3, simbetas3, ci=0.95) %>% as.data.frame() %>% mutate(iv='Economic Deprivation', dv='Populist Voting', Ideology='Left')
out4 <- logitsimfd(mm4, simbetas4, ci=0.95) %>% as.data.frame() %>% mutate(iv='Power Deprivation', dv='Populist Voting', Ideology='Left')
out5 <- logitsimfd(mm5, simbetas5, ci=0.95) %>% as.data.frame() %>% mutate(iv='Deprivation Scale', dv='Populist Voting', Ideology='Right')
out6 <- logitsimfd(mm6, simbetas6, ci=0.95) %>% as.data.frame() %>% mutate(iv='Social Deprivation', dv='Populist Voting', Ideology='Right')
out7 <- logitsimfd(mm7, simbetas7, ci=0.95) %>% as.data.frame() %>% mutate(iv='Economic Deprivation', dv='Populist Voting', Ideology='Right')
out8 <- logitsimfd(mm8, simbetas8, ci=0.95) %>% as.data.frame() %>% mutate(iv='Power Deprivation', dv='Populist Voting', Ideology='Right')

populist_voting_out <- bind_rows(out1, out2, out3, out4, out5, out6, out7, out8) 

plot2 <- bind_rows(populist_attitudes_out, populist_voting_out) %>%
  mutate(iv = factor(iv, levels=rev(c('Deprivation Scale','Social Deprivation','Economic Deprivation', 'Power Deprivation'))),
         Ideology=factor(Ideology)) %>%
  ggplot(aes(x=iv, y=pe, ymin=lower, ymax=upper, color=Ideology, shape=Ideology)) +
  coord_flip() + 
  geom_hline(yintercept=0, linetype=2, color='red') +
  geom_errorbar(width=0.1, position=position_dodge(width=0.5)) +
  geom_point(size=7, fill='white', position=position_dodge(width=0.5)) +
  facet_wrap(~dv) +
  labs(x='',y='Change in Y\n(min-max Deprivation)') +
  theme_bw() +
  geom_text(aes(x=iv, y=pe, label=round(pe,2)),size=2, position=position_dodge(width=0.5),show.legend = FALSE) +
  scale_color_manual(values=c('Blue','Red')) +
  scale_shape_manual(values=c(21,22))
plot2
ggsave(plot2, width=6, height=3.5, 
       file='figs/figure3.pdf')

## MODELS: by East vs. West ----

demo1e <- "total_dep_s*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ cze + est + hun + lva + ltu + pol + rou + svk + svn"
demo2e <- "social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ cze + est + hun + lva + ltu + pol + rou + svk + svn"
demo3e <- "econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ cze + est + hun + lva + ltu + pol + rou + svk + svn"
demo4e <- "power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ cze + est + hun + lva + ltu + pol + rou + svk + svn"

demo1w <- "total_dep_s*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ dnk + fra + deu + ita + nld + esp + swe"
demo2w <- "social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ dnk + fra + deu + ita + nld + esp + swe"
demo3w <- "econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ dnk + fra + deu + ita + nld + esp + swe"
demo4w <- "power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ dnk + fra + deu + ita + nld + esp + swe"

east <- df_t %>% filter(east==1) %>% as_survey_design(ids=country,weights=weighttot)
west <- df_t %>% filter(west==1) %>% as_survey_design(ids=country,weights=weighttot)

model1 <- svyglm(formula(paste0("populism_scale ~ ",demo1e)), east, family=gaussian(link='identity'))
model2 <- svyglm(formula(paste0("populism_scale ~ ",demo2e)), east, family=gaussian(link='identity'))
model3 <- svyglm(formula(paste0("populism_scale ~ ",demo3e)), east, family=gaussian(link='identity'))
model4 <- svyglm(formula(paste0("populism_scale ~ ",demo4e)), east, family=gaussian(link='identity'))
model5 <- svyglm(formula(paste0("populism_scale ~ ",demo1w)), west, family=gaussian(link='identity'))
model6 <- svyglm(formula(paste0("populism_scale ~ ",demo2w)), west, family=gaussian(link='identity'))
model7 <- svyglm(formula(paste0("populism_scale ~ ",demo3w)), west, family=gaussian(link='identity'))
model8 <- svyglm(formula(paste0("populism_scale ~ ",demo4w)), west, family=gaussian(link='identity'))
texreg::texreg(list(model1, model2, model3, model4, model5, model6, model7, model8),                   
               omit.coef="(deu)|(ita)|(esp)|(nld)|(fra)|(swe)|(rou)|(hun)|(pol)|(est)|(aut)|(dnk)|(cze)|(lva)|(ltu)|(svn)|(svk)|(bgr)",
               naive=T,reorder.coef = c(1,2,21,23,25,3:19,20,22,24,26),
               custom.coef.names = c('Intercept','Deprivation Scale','Incumbent','Female','Professional',
                                     'Service Worker','Worker','Age 35-49','Age 50-59','Age 60 plus',
                                     'Med Income','High Income','Med Education','High Education','Married',
                                     'Foreign Born','Worried COVID','Unemployed','Religious','Dep Scale * Incumbent',
                                     'Social Deprivation','Soc Dep * Incumbent','Economic Deprivation', 'Econ Dep * Incumbent',
                                     'Power Deprivation','Power Dep * Incumbent'),custom.note = '',
               custom.header = list('Populist Attitudes (East)'=1:4,'Populist Attitudes (West)'=5:8),
               file = 'tabs/tableB9.tex')
set.seed(12345)
pe <- coef(model1); vc <- vcov(model1)
simbetas1 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model2); vc <- vcov(model2)
simbetas2 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model3); vc <- vcov(model3)
simbetas3 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model4); vc <- vcov(model4)
simbetas4 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model5); vc <- vcov(model5)
simbetas5 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model6); vc <- vcov(model6)
simbetas6 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model7); vc <- vcov(model7)
simbetas7 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model8); vc <- vcov(model8)
simbetas8 <- MASS::mvrnorm(10000, pe, vc)

mm1 <- cfMake(populism_scale ~ total_dep_s*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + cze + est + hun + lva + ltu + pol + rou + svk + svn, data = df_t[df_t$east==1,], nscen = 1)
mm2 <- cfMake(populism_scale ~ social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + cze + est + hun + lva + ltu + pol + rou + svk + svn, data = df_t[df_t$east==1,], nscen = 1)
mm3 <- cfMake(populism_scale ~ econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + cze + est + hun + lva + ltu + pol + rou + svk + svn, data = df_t[df_t$east==1,], nscen = 1)
mm4 <- cfMake(populism_scale ~ power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + cze + est + hun + lva + ltu + pol + rou + svk + svn, data = df_t[df_t$east==1,], nscen = 1)
mm5 <- cfMake(populism_scale ~ total_dep_s*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ dnk + fra + deu + ita + nld + esp + swe, data = df_t[df_t$west==1,], nscen = 1)
mm6 <- cfMake(populism_scale ~ social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ dnk + fra + deu + ita + nld + esp + swe, data = df_t[df_t$west==1,], nscen = 1)
mm7 <- cfMake(populism_scale ~ econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ dnk + fra + deu + ita + nld + esp + swe, data = df_t[df_t$west==1,], nscen = 1)
mm8 <- cfMake(populism_scale ~ power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ dnk + fra + deu + ita + nld + esp + swe, data = df_t[df_t$west==1,], nscen = 1)

mm1 <- cfChange(mm1, 'total_dep_s',xpre=0, x=1, scen=1)
mm2 <- cfChange(mm2, 'social_dep',xpre=0, x=1, scen=1)
mm3 <- cfChange(mm3, 'econ_dep',xpre=0, x=1, scen=1)
mm4 <- cfChange(mm4, 'power_dep',xpre=0, x=1, scen=1)
mm5 <- cfChange(mm5, 'total_dep_s',xpre=0, x=1, scen=1)
mm6 <- cfChange(mm6, 'social_dep',xpre=0, x=1, scen=1)
mm7 <- cfChange(mm7, 'econ_dep',xpre=0, x=1, scen=1)
mm8 <- cfChange(mm8, 'power_dep',xpre=0, x=1, scen=1)

mm1 <- cfChange(mm1, 'incumbent',xpre=0, x=0, scen=1)
mm2 <- cfChange(mm2, 'incumbent',xpre=0, x=0, scen=1)
mm3 <- cfChange(mm3, 'incumbent',xpre=0, x=0, scen=1)
mm4 <- cfChange(mm4, 'incumbent',xpre=0, x=0, scen=1)
mm5 <- cfChange(mm5, 'incumbent',xpre=0, x=0, scen=1)
mm6 <- cfChange(mm6, 'incumbent',xpre=0, x=0, scen=1)
mm7 <- cfChange(mm7, 'incumbent',xpre=0, x=0, scen=1)
mm8 <- cfChange(mm8, 'incumbent',xpre=0, x=0, scen=1)

out1 <- linearsimfd(mm1, simbetas1, ci=0.95) %>% as.data.frame() %>% mutate(iv='Deprivation Scale', dv='Populist Attitudes', Region='East')
out2 <- linearsimfd(mm2, simbetas2, ci=0.95) %>% as.data.frame() %>% mutate(iv='Social Deprivation', dv='Populist Attitudes', Region='East')
out3 <- linearsimfd(mm3, simbetas3, ci=0.95) %>% as.data.frame() %>% mutate(iv='Economic Deprivation', dv='Populist Attitudes', Region='East')
out4 <- linearsimfd(mm4, simbetas4, ci=0.95) %>% as.data.frame() %>% mutate(iv='Power Deprivation', dv='Populist Attitudes', Region='East')
out5 <- linearsimfd(mm5, simbetas5, ci=0.95) %>% as.data.frame() %>% mutate(iv='Deprivation Scale', dv='Populist Attitudes', Region='West')
out6 <- linearsimfd(mm6, simbetas6, ci=0.95) %>% as.data.frame() %>% mutate(iv='Social Deprivation', dv='Populist Attitudes', Region='West')
out7 <- linearsimfd(mm7, simbetas7, ci=0.95) %>% as.data.frame() %>% mutate(iv='Economic Deprivation', dv='Populist Attitudes', Region='West')
out8 <- linearsimfd(mm8, simbetas8, ci=0.95) %>% as.data.frame() %>% mutate(iv='Power Deprivation', dv='Populist Attitudes', Region='West')

populist_attitudes_out <- bind_rows(out1, out2, out3, out4, out5, out6, out7, out8) 

model1 <- svyglm(formula(paste0("voted_populist ~ ",demo1e)), east, family=binomial(link='logit'))
model2 <- svyglm(formula(paste0("voted_populist ~ ",demo2e)), east, family=binomial(link='logit'))
model3 <- svyglm(formula(paste0("voted_populist ~ ",demo3e)), east, family=binomial(link='logit'))
model4 <- svyglm(formula(paste0("voted_populist ~ ",demo4e)), east, family=binomial(link='logit'))
model5 <- svyglm(formula(paste0("voted_populist ~ ",demo1w)), west, family=binomial(link='logit'))
model6 <- svyglm(formula(paste0("voted_populist ~ ",demo2w)), west, family=binomial(link='logit'))
model7 <- svyglm(formula(paste0("voted_populist ~ ",demo3w)), west, family=binomial(link='logit'))
model8 <- svyglm(formula(paste0("voted_populist ~ ",demo4w)), west, family=binomial(link='logit'))
texreg::texreg(list(model1, model2, model3, model4, model5, model6, model7, model8),                   
               omit.coef="(deu)|(ita)|(esp)|(nld)|(fra)|(swe)|(rou)|(hun)|(pol)|(est)|(aut)|(dnk)|(cze)|(lva)|(ltu)|(svn)|(svk)|(bgr)",
               naive=T,reorder.coef = c(1,2,21,23,25,3:19,20,22,24,26),
               custom.coef.names = c('Intercept','Deprivation Scale','Incumbent','Female','Professional',
                                     'Service Worker','Worker','Age 35-49','Age 50-59','Age 60 plus',
                                     'Med Income','High Income','Med Education','High Education','Married',
                                     'Foreign Born','Worried COVID','Unemployed','Religious','Dep Scale * Incumbent',
                                     'Social Deprivation','Soc Dep * Incumbent','Economic Deprivation', 'Econ Dep * Incumbent',
                                     'Power Deprivation','Power Dep * Incumbent'),custom.note = '',
               custom.header = list('Populist Voting (East)'=1:4,'Populist Voting (West)'=5:8),
               file = 'tabs/tableB10.tex')
set.seed(12345)
pe <- coef(model1); vc <- vcov(model1)
simbetas1 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model2); vc <- vcov(model2)
simbetas2 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model3); vc <- vcov(model3)
simbetas3 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model4); vc <- vcov(model4)
simbetas4 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model5); vc <- vcov(model5)
simbetas5 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model6); vc <- vcov(model6)
simbetas6 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model7); vc <- vcov(model7)
simbetas7 <- MASS::mvrnorm(10000, pe, vc)
pe <- coef(model8); vc <- vcov(model8)
simbetas8 <- MASS::mvrnorm(10000, pe, vc)

mm1 <- cfMake(voted_populist ~ total_dep_s*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + cze + est + hun + lva + ltu + pol + rou + svk + svn, data = df_t[df_t$east==1,], nscen = 1)
mm2 <- cfMake(voted_populist ~ social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + cze + est + hun + lva + ltu + pol + rou + svk + svn, data = df_t[df_t$east==1,], nscen = 1)
mm3 <- cfMake(voted_populist ~ econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + cze + est + hun + lva + ltu + pol + rou + svk + svn, data = df_t[df_t$east==1,], nscen = 1)
mm4 <- cfMake(voted_populist ~ power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + cze + est + hun + lva + ltu + pol + rou + svk + svn, data = df_t[df_t$east==1,], nscen = 1)
mm5 <- cfMake(voted_populist ~ total_dep_s*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ dnk + fra + deu + ita + nld + esp + swe, data = df_t[df_t$west==1,], nscen = 1)
mm6 <- cfMake(voted_populist ~ social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ dnk + fra + deu + ita + nld + esp + swe, data = df_t[df_t$west==1,], nscen = 1)
mm7 <- cfMake(voted_populist ~ econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ dnk + fra + deu + ita + nld + esp + swe, data = df_t[df_t$west==1,], nscen = 1)
mm8 <- cfMake(voted_populist ~ power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ dnk + fra + deu + ita + nld + esp + swe, data = df_t[df_t$west==1,], nscen = 1)

mm1 <- cfChange(mm1, 'total_dep_s',xpre=0, x=1, scen=1)
mm2 <- cfChange(mm2, 'social_dep',xpre=0, x=1, scen=1)
mm3 <- cfChange(mm3, 'econ_dep',xpre=0, x=1, scen=1)
mm4 <- cfChange(mm4, 'power_dep',xpre=0, x=1, scen=1)
mm5 <- cfChange(mm5, 'total_dep_s',xpre=0, x=1, scen=1)
mm6 <- cfChange(mm6, 'social_dep',xpre=0, x=1, scen=1)
mm7 <- cfChange(mm7, 'econ_dep',xpre=0, x=1, scen=1)
mm8 <- cfChange(mm8, 'power_dep',xpre=0, x=1, scen=1)

mm1 <- cfChange(mm1, 'incumbent',xpre=0, x=0, scen=1)
mm2 <- cfChange(mm2, 'incumbent',xpre=0, x=0, scen=1)
mm3 <- cfChange(mm3, 'incumbent',xpre=0, x=0, scen=1)
mm4 <- cfChange(mm4, 'incumbent',xpre=0, x=0, scen=1)
mm5 <- cfChange(mm5, 'incumbent',xpre=0, x=0, scen=1)
mm6 <- cfChange(mm6, 'incumbent',xpre=0, x=0, scen=1)
mm7 <- cfChange(mm7, 'incumbent',xpre=0, x=0, scen=1)
mm8 <- cfChange(mm8, 'incumbent',xpre=0, x=0, scen=1)

out1 <- logitsimfd(mm1, simbetas1, ci=0.95) %>% as.data.frame() %>% mutate(iv='Deprivation Scale', dv='Populist Voting', Region='East')
out2 <- logitsimfd(mm2, simbetas2, ci=0.95) %>% as.data.frame() %>% mutate(iv='Social Deprivation', dv='Populist Voting', Region='East')
out3 <- logitsimfd(mm3, simbetas3, ci=0.95) %>% as.data.frame() %>% mutate(iv='Economic Deprivation', dv='Populist Voting', Region='East')
out4 <- logitsimfd(mm4, simbetas4, ci=0.95) %>% as.data.frame() %>% mutate(iv='Power Deprivation', dv='Populist Voting', Region='East')
out5 <- logitsimfd(mm5, simbetas5, ci=0.95) %>% as.data.frame() %>% mutate(iv='Deprivation Scale', dv='Populist Voting', Region='West')
out6 <- logitsimfd(mm6, simbetas6, ci=0.95) %>% as.data.frame() %>% mutate(iv='Social Deprivation', dv='Populist Voting', Region='West')
out7 <- logitsimfd(mm7, simbetas7, ci=0.95) %>% as.data.frame() %>% mutate(iv='Economic Deprivation', dv='Populist Voting', Region='West')
out8 <- logitsimfd(mm8, simbetas8, ci=0.95) %>% as.data.frame() %>% mutate(iv='Power Deprivation', dv='Populist Voting', Region='West')

populist_voting_out <- bind_rows(out1, out2, out3, out4, out5, out6, out7, out8) 

plot3 <- bind_rows(populist_attitudes_out, populist_voting_out) %>%
  mutate(iv = factor(iv, levels=rev(c('Deprivation Scale','Social Deprivation','Economic Deprivation', 'Power Deprivation'))),
         Region=factor(Region)) %>%
  ggplot(aes(x=iv, y=pe, ymin=lower, ymax=upper, color=Region, shape=Region)) +
  coord_flip() + 
  geom_hline(yintercept=0, linetype=2, color='red') +
  geom_errorbar(width=0.1, position=position_dodge(width=0.5)) +
  geom_point(size=7, fill='white', position=position_dodge(width=0.5)) +
  facet_wrap(~dv) +
  labs(x='',y='Change in Y\n(min-max Deprivation)') +
  theme_bw() +
  geom_text(aes(x=iv, y=pe, label=round(pe,2)),size=2, position=position_dodge(width=0.5),show.legend = FALSE) +
  scale_color_manual(values=c('Blue','Red')) +
  scale_shape_manual(values=c(21,22))
plot3
ggsave(plot3, width=6, height=3.5, 
       file='figs/figure4_top.pdf')

# robustness check countries with robust left and right wing parties
countries <- df_t %>%
  group_by(country, leftist) %>%
  summarise(y = weighted.mean(voted_populist, na.rm=T, w=weight)) %>%
  filter(y > 0.15) %>%
  as.data.frame() %>%
  filter(!country %in% c(
    'Austria (AUT)','Hungary (HUN)','Poland (POL)'
  )) %>%
  ungroup() %>%
  select(country) %>%
  unique() %>% pull()

demo1 <- "total_dep_s*incumbent + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious + deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"
demo2 <- "social_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"
demo3 <- "econ_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"
demo4 <- "power_dep*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"

my_des <- df_t  %>% filter(country %in% countries) %>% as_survey_design(ids=country,weights=weighttot)

model1 <- svyglm(formula(paste0("populism_scale ~ ",demo1)), my_des, family=gaussian(link='identity'))
model2 <- svyglm(formula(paste0("populism_scale ~ ",demo2)), my_des, family=gaussian(link='identity'))
model3 <- svyglm(formula(paste0("populism_scale ~ ",demo3)), my_des, family=gaussian(link='identity'))
model4 <- svyglm(formula(paste0("populism_scale ~ ",demo4)), my_des, family=gaussian(link='identity'))
model5 <- svyglm(formula(paste0("voted_populist ~ ",demo1)), my_des, family=binomial(link='logit'))
model6 <- svyglm(formula(paste0("voted_populist ~ ",demo2)), my_des, family=binomial(link='logit'))
model7 <- svyglm(formula(paste0("voted_populist ~ ",demo3)), my_des, family=binomial(link='logit'))
model8 <- svyglm(formula(paste0("voted_populist ~ ",demo4)), my_des, family=binomial(link='logit'))

texreg::texreg(list(model1, model2, model3, model4, model5, model6, model7, model8),                   
               omit.coef="(deu)|(ita)|(esp)|(nld)|(fra)|(swe)|(rou)|(hun)|(pol)|(est)|(aut)|(dnk)|(cze)|(lva)|(ltu)|(svn)|(svk)|(bgr)",
               naive=T,reorder.coef = c(1,2,21,23,25,3:19,20,22,24,26),
               custom.coef.names = c('Intercept','Deprivation Scale','Incumbent','Female','Professional',
                                     'Service Worker','Worker','Age 35-49','Age 50-59','Age 60 plus',
                                     'Med Income','High Income','Med Education','High Education','Married',
                                     'Foreign Born','Worried COVID','Unemployed','Religious','Dep Scale * Incumbent',
                                     'Social Deprivation','Soc Dep * Incumbent','Economic Deprivation', 'Econ Dep * Incumbent',
                                     'Power Deprivation','Power Dep * Incumbent'),custom.note = '',
               custom.header = list('Populist Attitudes'=1:4,'Populist Voting'=5:8),
               file = 'tabs/tableB11.tex')


# Results country by country

# list of all countries
countries <- unique(df_t$country)

# formulas for incumbent and fb
demo1 <- "total_dep_s_new + incumbent + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious"
demo2 <- "social_dep_new + incumbent + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious"
demo3 <- "econ_dep_new + incumbent + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious"
demo4 <- "power_dep_new + incumbent + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious"

# formulas for incumbent no foreign
demo1_nf <- "total_dep_s_new + incumbent + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + worried_corona + unemployed + religious"
demo2_nf <- "social_dep_new + incumbent + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + worried_corona + unemployed + religious"
demo3_nf <- "econ_dep_new + incumbent + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + worried_corona + unemployed + religious"
demo4_nf <- "power_dep_new + incumbent + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + worried_corona + unemployed + religious"

# formulas for no incum but foreign
demo1_i_nf  <- "total_dep_s_new + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious"
demo2_i_nf  <- "social_dep_new + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious"
demo3_i_nf  <- "econ_dep_new + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious"
demo4_i_nf  <- "power_dep_new + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious"

# formulas no incumbent no foreign
demo1_ni_nf <- "total_dep_s_new + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + worried_corona + unemployed + religious"
demo2_ni_nf <- "social_dep_new + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + worried_corona + unemployed + religious"
demo3_ni_nf <- "econ_dep_new + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + worried_corona + unemployed + religious"
demo4_ni_nf <- "power_dep_new + female + professionals + service_workers + worker + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + worried_corona + unemployed + religious"

extractCoefs <- function(model1, country, outcome){
  broom::tidy(model1) %>%
    filter(term %in% c('total_dep_s_new','social_dep_new','econ_dep_new','power_dep_new')) %>%
    mutate(country = country,
           l = estimate - 1.96 * std.error,
           u = estimate + 1.96 * std.error,
           outcome=outcome) 
}

extractFDs <- function(country_input){
  
  # incumbent 
  incumbent = ifelse(sum(df_t$incumbent[df_t$country == country_input])>10,1,0)
  fb = ifelse(mean(df_t$foreign_born[df_t$country == country_input])>0.02,1,0)
  
  if(incumbent == 1 & fb == 1){
    
    form1 <- formula(paste0("voted_populist ~ ",demo1))
    form2 <- formula(paste0("voted_populist ~ ",demo2))
    form3 <- formula(paste0("voted_populist ~ ",demo3))
    form4 <- formula(paste0("voted_populist ~ ",demo4))
    
  } else if(incumbent == 1 & fb == 0) {
    
    form1 <- formula(paste0("voted_populist ~ ",demo1_nf))
    form2 <- formula(paste0("voted_populist ~ ",demo2_nf))
    form3 <- formula(paste0("voted_populist ~ ",demo3_nf))
    form4 <- formula(paste0("voted_populist ~ ",demo4_nf))
    
  } else if(incumbent == 0 & fb == 0){
    
    form1 <- formula(paste0("voted_populist ~ ",demo1_ni_nf))
    form2 <- formula(paste0("voted_populist ~ ",demo2_ni_nf))
    form3 <- formula(paste0("voted_populist ~ ",demo3_ni_nf))
    form4 <- formula(paste0("voted_populist ~ ",demo4_ni_nf))
    
  } else if(incumbent == 0 & fb == 1){
    
    form1 <- formula(paste0("voted_populist ~ ",demo1_i_nf))
    form2 <- formula(paste0("voted_populist ~ ",demo2_i_nf))
    form3 <- formula(paste0("voted_populist ~ ",demo3_i_nf))
    form4 <- formula(paste0("voted_populist ~ ",demo4_i_nf))
    
  }
  
  my_des <- df_t  %>% filter(country %in% country_input)
  
  res.out1 <- glm(form1, data=my_des, family=binomial(link='logit'))
  res.out2 <- glm(form2, data=my_des, family=binomial(link='logit'))
  res.out3 <- glm(form3, data=my_des, family=binomial(link='logit'))
  res.out4 <- glm(form4, data=my_des, family=binomial(link='logit'))
  
  summary(res.out1)
  summary(res.out2)
  summary(res.out3)
  summary(res.out4)
  
  mmout1 <- cfMake(form1, my_des, nscen=1)
  mmout2 <- cfMake(form2, my_des, nscen=1)
  mmout3 <- cfMake(form3, my_des, nscen=1)
  mmout4 <- cfMake(form4, my_des, nscen=1)
  
  if(incumbent == 1 & fb == 1){
    mmout1 <- cfChange(mmout1, 'total_dep_s_new',xpre=0, x=1, scen=1)
    mmout2 <- cfChange(mmout2, 'social_dep_new',xpre=0, x=1, scen=1)
    mmout3 <- cfChange(mmout3, 'econ_dep_new',xpre=0, x=1, scen=1)
    mmout4 <- cfChange(mmout4, 'power_dep_new',xpre=0, x=1, scen=1)
    
    mmout1 <- cfChange(mmout1, 'incumbent',xpre=0, x=0, scen=1)
    mmout2 <- cfChange(mmout2, 'incumbent',xpre=0, x=0, scen=1)
    mmout3 <- cfChange(mmout3, 'incumbent',xpre=0, x=0, scen=1)
    mmout4 <- cfChange(mmout4, 'incumbent',xpre=0, x=0, scen=1)
    
  } else if(incumbent == 1 & fb == 0) {
    mmout1 <- cfChange(mmout1, 'total_dep_s_new',xpre=0, x=1, scen=1)
    mmout2 <- cfChange(mmout2, 'social_dep_new',xpre=0, x=1, scen=1)
    mmout3 <- cfChange(mmout3, 'econ_dep_new',xpre=0, x=1, scen=1)
    mmout4 <- cfChange(mmout4, 'power_dep_new',xpre=0, x=1, scen=1)
    
    mmout1 <- cfChange(mmout1, 'incumbent',xpre=0, x=0, scen=1)
    mmout2 <- cfChange(mmout2, 'incumbent',xpre=0, x=0, scen=1)
    mmout3 <- cfChange(mmout3, 'incumbent',xpre=0, x=0, scen=1)
    mmout4 <- cfChange(mmout4, 'incumbent',xpre=0, x=0, scen=1)  
  } else if(incumbent == 0 & fb == 0){
    mmout1 <- cfChange(mmout1, 'total_dep_s_new',xpre=0, x=1, scen=1)
    mmout2 <- cfChange(mmout2, 'social_dep_new',xpre=0, x=1, scen=1)
    mmout3 <- cfChange(mmout3, 'econ_dep_new',xpre=0, x=1, scen=1)
    mmout4 <- cfChange(mmout4, 'power_dep_new',xpre=0, x=1, scen=1)
  } else if(incumbent == 0 & fb == 1){
    mmout1 <- cfChange(mmout1, 'total_dep_s_new',xpre=0, x=1, scen=1)
    mmout2 <- cfChange(mmout2, 'social_dep_new',xpre=0, x=1, scen=1)
    mmout3 <- cfChange(mmout3, 'econ_dep_new',xpre=0, x=1, scen=1)
    mmout4 <- cfChange(mmout4, 'power_dep_new',xpre=0, x=1, scen=1)
  }
  
  out1 = marginaleffects::avg_comparisons(res.out1,
                                          newdata = mmout1$xpre,
                                          variables = list(total_dep_s_new = c(0,1)))
  out2 = marginaleffects::avg_comparisons(res.out2,
                                          newdata = mmout2$xpre,
                                          variables = list(social_dep_new = c(0,1)))
  out3 = marginaleffects::avg_comparisons(res.out3,
                                          newdata = mmout3$xpre,
                                          variables = list(econ_dep_new = c(0,1)))
  out4 = marginaleffects::avg_comparisons(res.out4,
                                          newdata = mmout4$xpre,
                                          variables = list(power_dep_new = c(0,1)))
  
  res.all <- bind_rows(out1, out2, out3, out4) %>%
    as_tibble() %>%
    mutate(country = country_input,
           outcome='Populist Voting') 
  
  return(res.all)
}

results_inc = vector('list', length(countries))

for(i in 1:length(countries)){
  
  incumbent = ifelse(sum(df_t$incumbent[df_t$country == countries[i]])>10,1,0)
  fb = ifelse(sum(df_t$foreign_born[df_t$country == countries[i]])>10,1,0)
  
  my_des <- df_t  %>% 
    filter(country %in% countries[i])
  
  if(incumbent == 1 & fb == 1){
    
    model1 <- lm_robust(formula(paste0("populism_scale ~ ",demo1)), my_des)
    model2 <- lm_robust(formula(paste0("populism_scale ~ ",demo2)), my_des)
    model3 <- lm_robust(formula(paste0("populism_scale ~ ",demo3)), my_des)
    model4 <- lm_robust(formula(paste0("populism_scale ~ ",demo4)), my_des)
    
  } else if(incumbent == 1 & fb == 0) {
    
    model1 <- lm_robust(formula(paste0("populism_scale ~ ",demo1_nf)), my_des)
    model2 <- lm_robust(formula(paste0("populism_scale ~ ",demo2_nf)), my_des)
    model3 <- lm_robust(formula(paste0("populism_scale ~ ",demo3_nf)), my_des)
    model4 <- lm_robust(formula(paste0("populism_scale ~ ",demo4_nf)), my_des)
    
  } else if(incumbent == 0 & fb == 0){
    
    model1 <- lm_robust(formula(paste0("populism_scale ~ ",demo1_ni_nf)), my_des)
    model2 <- lm_robust(formula(paste0("populism_scale ~ ",demo2_ni_nf)), my_des)
    model3 <- lm_robust(formula(paste0("populism_scale ~ ",demo3_ni_nf)), my_des)
    model4 <- lm_robust(formula(paste0("populism_scale ~ ",demo4_ni_nf)), my_des)
    
  } else if(incumbent == 0 & fb == 1){
    
    model1 <- lm_robust(formula(paste0("populism_scale ~ ",demo1_i_nf)), my_des)
    model2 <- lm_robust(formula(paste0("populism_scale ~ ",demo2_i_nf)), my_des)
    model3 <- lm_robust(formula(paste0("populism_scale ~ ",demo3_i_nf)), my_des)
    model4 <- lm_robust(formula(paste0("populism_scale ~ ",demo4_i_nf)), my_des)
    
  }
  
  results_attitudes <- bind_rows(
    extractCoefs(model1, countries[i], 'Populist Attitudes'),
    extractCoefs(model2, countries[i], 'Populist Attitudes'),
    extractCoefs(model3, countries[i], 'Populist Attitudes'),
    extractCoefs(model4, countries[i], 'Populist Attitudes')
  )
  
  results_voting <- extractFDs(countries[i]) %>%
    rename(l = conf.low, u=conf.high) %>%
    mutate(outcome='Populist Voting')
  
  results_inc[[i]] <- bind_rows(results_attitudes, results_voting) %>%
    select(term, estimate, l, u, country, outcome)
  
  print(i)
}

results_inc <- bind_rows(results_inc) %>%
  mutate(outcome = case_when(
    outcome == 'populism_scale' ~'Populist Attitudes',
    TRUE ~ outcome
  )) 
p1 <- results_inc %>%
  filter(outcome=='Populist Attitudes') %>%
  mutate(term = case_when(
    term == 'total_dep_s_new' ~ 'Total Deprivation',
    term == 'social_dep_new' ~ 'Social Deprivation',
    term == 'power_dep_new' ~ 'Power Deprivation',
    term == 'econ_dep_new' ~ 'Economic Deprivation'
  )) %>%
  ggplot(aes(x=term, y=estimate, ymin=l, ymax=u)) +
  geom_hline(yintercept=0, color='red', linetype=2) +
  geom_errorbar(width=0.1) +
  geom_point(shape=21, fill='white') +
  facet_wrap(~country) +
  coord_flip() +
  theme_bw()+
  labs(x='', y='Change in Populist Attitudes')
p1

# 
p2 <- results_inc %>%
  filter(outcome=='Populist Voting') %>%
  filter(!country %in% c('Romania (ROU)')) %>%
  mutate(term = case_when(
    term == 'total_dep_s_new' ~ 'Total Deprivation',
    term == 'social_dep_new' ~ 'Social Deprivation',
    term == 'power_dep_new' ~ 'Power Deprivation',
    term == 'econ_dep_new' ~ 'Economic Deprivation')) %>%
  ggplot(aes(x=term, y=estimate, ymin=l, ymax=u)) +
  geom_hline(yintercept=0, color='red', linetype=2) +
  geom_errorbar(width=0.1) +
  geom_point(shape=21, fill='white') +
  facet_wrap(~country) +
  coord_flip() +
  theme_bw() +
  labs(x='', y='Change in Pr(Vote Populist)')
p2

p1 / p2
ggsave(width=10, height=9,
       filename = 'figs/figureC3.pdf')


# Figure 4 bottom

bothdat <- results_inc %>%
  mutate(term = case_when(
    term == 'total_dep_s_new' ~ 'Total Deprivation',
    term == 'social_dep_new' ~ 'Social Deprivation',
    term == 'power_dep_new' ~ 'Power Deprivation',
    term == 'econ_dep_new' ~ 'Economic Deprivation'
  )) %>%
  filter(term == 'Total Deprivation') %>%
  mutate(region = ifelse(
    country %in% c('Bulgaria (BGR)','Czech Republic (CZE)',
                   'Estonia (EST)','Hungary (HUN)','Latvia (LVA)',
                   'Lithuania (LTU)','Poland (POL)','Romania (ROU)',
                   'Slovakia (SVK)','Slovenia (SVN)'),'East','West'
  ), region=factor(region, levels=c('West','East')))
top <- bothdat %>%
  filter(region == 'West') %>%
  ggplot(aes(x=reorder(country,estimate), y=estimate, ymin=l, ymax=u)) +
  geom_hline(yintercept=0, color='red', linetype=2) +
  geom_pointrange(position=position_dodge(width=0.5), fill='white',shape=21) +
  facet_wrap(~outcome) +
  coord_flip() +
  theme_bw()+
  labs(x='', y='',title = 'West') +
  scale_shape_manual(values = c(21,24),name='') +
  scale_color_brewer(palette='Set2',name='') +
  scale_y_continuous(limits=c(-0.75,1.1)) +
  theme(plot.title = element_text(hjust = 0.5))

bottom_dat <- bothdat %>%
  filter(region == 'East') 
bottom <- bottom_dat %>%
  ggplot(aes(x=reorder(country,estimate), y=estimate, ymin=l, ymax=u)) +
  geom_hline(yintercept=0, color='red', linetype=2) +
  geom_pointrange(position=position_dodge(width=0.5), fill='white',shape=21) +
  facet_wrap(~outcome) +
  coord_flip() +
  theme_bw()+
  labs(x='', y='',title = 'East') +
  scale_shape_manual(values = c(21,24),name='') +
  scale_color_brewer(palette='Set2',name='') +
  scale_y_continuous(limits=c(-0.75,1.1)) +
  theme(plot.title = element_text(hjust = 0.5))

top / bottom
ggsave(width=8, height=6,
       filename = 'figs/figure4_bottom.pdf')


# alternate specifications ----
df_t$social_dep_factor <- NA
df_t$social_dep_factor[df_t$social_dep > 0.5] <- 'Negative'
df_t$social_dep_factor[df_t$social_dep == 0.5] <- 'Neutral'
df_t$social_dep_factor[df_t$social_dep < 0.5] <- 'Positive'
df_t$social_dep_factor <- factor(df_t$social_dep_factor, levels=c(
  'Neutral','Negative','Positive'
))
table(df_t$social_dep_factor)

df_t$econ_dep_factor <- NA
df_t$econ_dep_factor[df_t$econ_dep > 0.5] <- 'Negative'
df_t$econ_dep_factor[df_t$econ_dep == 0.5] <- 'Neutral'
df_t$econ_dep_factor[df_t$econ_dep < 0.5] <- 'Positive'
df_t$econ_dep_factor <- factor(df_t$econ_dep_factor, levels=c(
  'Neutral','Negative','Positive'
))
table(df_t$econ_dep_factor)

df_t$power_dep_factor <- NA
df_t$power_dep_factor[df_t$power_dep > 0.5] <- 'Negative'
df_t$power_dep_factor[df_t$power_dep == 0.5] <- 'Neutral'
df_t$power_dep_factor[df_t$power_dep < 0.5] <- 'Positive'
df_t$power_dep_factor <- factor(df_t$power_dep_factor, levels=c(
  'Neutral','Negative','Positive'
))

my_des <- df_t  %>% as_survey_design(ids=country,weights=weighttot)

demo1 <- "social_dep_factor*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"
demo2 <- "econ_dep_factor*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"
demo3 <- "power_dep_factor*incumbent + female + professionals + service_workers + worker  + age_35_49 + age_50_59 + age_60 + inc_med + inc_high + edu_med + edu_high + married + foreign_born + worried_corona + unemployed + religious+ deu + ita + esp + nld + fra + swe + rou + hun + pol + est + aut + dnk + cze + lva + ltu + svn + svk + bgr"

model1 <- svyglm(formula(paste0("populism_scale ~ ",demo1)), my_des, family=gaussian(link='identity'))
model2 <- svyglm(formula(paste0("populism_scale ~ ",demo2)), my_des, family=gaussian(link='identity'))
model3 <- svyglm(formula(paste0("populism_scale ~ ",demo3)), my_des, family=gaussian(link='identity'))
model4 <- svyglm(formula(paste0("voted_populist ~ ",demo1)), my_des, family=binomial(link='logit'))
model5 <- svyglm(formula(paste0("voted_populist ~ ",demo2)), my_des, family=binomial(link='logit'))
model6 <- svyglm(formula(paste0("voted_populist ~ ",demo3)), my_des, family=binomial(link='logit'))
texreg::texreg(list(model1, model2, model3, model4, model5, model6),
               omit.coef="(deu)|(ita)|(esp)|(nld)|(fra)|(swe)|(rou)|(hun)|(pol)|(est)|(aut)|(dnk)|(cze)|(lva)|(ltu)|(svn)|(svk)|(bgr)",
               naive=T,reorder.coef = c(1,2,3,23,24,27,28,4:22,25,26,29,30),
               custom.coef.names = c('Intercept','Social Dep (Neg)','Social Dep (Pos)',
                                     'Incumbent','Female','Professional',
                                     'Service Worker','Worker','Age 35-49','Age 50-59','Age 60 plus',
                                     'Med Income','High Income','Med Education','High Education','Married',
                                     'Foreign Born','Worried COVID','Unemployed','Religious',
                                     'Social Dep (Neg) * Incumbent','Social Dep (Pos) * Incumbent',
                                     'Econ Dep (Neg)','Econ Dep (Pos)',
                                     'Econ Dep (Neg * Incumbent)','Econ Dep (Pos) * Incumbent',
                                     'Power Dep (Neg)','Power Dep (Pos)',
                                     'Power Dep (Neg) * Incumbent','Power Dep (Pos) * Incumbent'),
               custom.note = '',
               custom.header = list('Populist Attitudes'=1:3,'Populist Voting'=4:6),
               file='tabs/tableB4.tex')



